/*
 * gaussian_integration.cpp
 *
 *  Created on: Jun 19, 2011
 *      Author: nzhao
 */
#include "share/gaussian_integration.h"
double * GAUSSIAN_POINT;
double * GAUSSIAN_WEIGHT;
void read_gaussian_quadrature_weights()
{
	double minX, maxX;
	int minOrder, maxOrder, n, k, n_idx;
	double ** gs_point, ** gs_weight;

	char Gauss_file[200];
	strcpy(Gauss_file, WORKSHOP);
	strcat(Gauss_file, DEBUG_PATH);
	strcat(Gauss_file, "GaussInt.dat");

	ifstream gs(Gauss_file, ios::in | ios::binary);

	gs.read((char *)&minX, sizeof(double));
	gs.read((char *)&maxX, sizeof(double));
	gs.read((char *)&minOrder, sizeof(int));
	gs.read((char *)&maxOrder, sizeof(int));

	cout << "minX = " << minX << "; maxX = " << maxX << endl;
	cout << "minOrder = " << minOrder << "; maxOrder = " << maxOrder<< endl;

	gs_point = new double * [maxOrder - minOrder + 1];
	gs_weight = new double * [maxOrder - minOrder + 1];
	for(n = minOrder; n <= maxOrder; n++)
	{
		n_idx = n - minOrder;
		gs_point[n_idx] = new double [n];
		gs.read((char *)gs_point[n_idx], n*sizeof(double));
	}
	for(n = minOrder; n <= maxOrder; n++)
	{
		n_idx = n - minOrder;
		gs_weight[n_idx] = new double [n];
		gs.read((char *)gs_weight[n_idx], n*sizeof(double));
	}
	gs.close();

	GAUSSIAN_POINT = new double [GAUSS_INT_ORDER];
	GAUSSIAN_WEIGHT= new double [GAUSS_INT_ORDER];
	for(k = 0; k < GAUSS_INT_ORDER; k++)
	{
		GAUSSIAN_POINT[k] = gs_point[GAUSS_INT_ORDER - minOrder][k];
		GAUSSIAN_WEIGHT[k] = gs_weight[GAUSS_INT_ORDER - minOrder][k];
	}

	for(n = minOrder; n <= maxOrder; n++)
	{
		n_idx = n - minOrder;
		delete [] gs_point[n_idx];
		delete [] gs_weight[n_idx];
	}
	delete [] gs_point;
	delete [] gs_weight;

}
